%% integral_loc
function [points, weights] = integral_loc(V, int_type)
[int_points_ref, int_weights_ref] = integral_ref(int_type);
if 2 == size(V,1) && strcmp(int_type, "1d")
    AB = 0.5*(V(2,:)-V(1,:));
    weights = norm(AB)*int_weights_ref;
    points = AB.*(int_points_ref+1)+V(1,:);
elseif 3 == size(V,1) && strcmp(int_type, "2d")
    Jacobi = V(2:3,:)-V(1,:);
    if 2 == size(V,2)
        J_det = det(Jacobi);
    else
        AB = V(2,:)-V(1,:);
        AC = V(3,:)-V(1,:);
        J_det = norm(cross(AB, AC));
    end
    weights = abs(J_det)*int_weights_ref;
    points = int_points_ref*Jacobi+V(1,:);
elseif 4 == size(V,1) && strcmp(int_type, "3d")
    Jacobi = V(2:4,:)-V(1,:);
    J_det = det(Jacobi);
    weights = abs(J_det)*int_weights_ref;
    points = int_points_ref*Jacobi+V(1,:);
else
    error("Invalid gauss-lengdre type.");
end
end

%% integral_ref
function [points, weights] = integral_ref(int_type)
switch int_type
    case "1d" % gauss-lengdre(4)
        points = [-0.8611363115940525752239465;
            -0.3399810435848562648026658;
            0.3399810435848562648026658;
            0.8611363115940525752239465];
        weights = [0.3478548451374538573730639, 0.6521451548625461426269361, 0.6521451548625461426269361, 0.3478548451374538573730639];
    case "2d" % Dunavant(7)
        points = [1/3, 1/3;
            0.470142064105115, 0.470142064105115;
            0.470142064105115, 0.059715871789770;
            0.059715871789770, 0.470142064105115;
            0.101286507323456, 0.101286507323456;
            0.101286507323456, 0.797426985353087;
            0.797426985353087, 0.101286507323456];
        weights = [0.1125, 0.066197076394253, 0.066197076394253, 0.066197076394253, 0.062969590272414, 0.062969590272414, 0.062969590272414];

        % points = [0.11270166537925829786104259255808, 0.1;
        %     0.11270166537925829786104259255808, 0.44364916731037085106947870372096;
        %     0.11270166537925829786104259255808, 0.78729833462074172434341789994505;
        %     0.5, 0.056350832689629148930521296279039;
        %     0.5, 0.25;
        %     0.5, 0.44364916731037085106947870372096;
        %     0.88729833462074170213895740744192, 0.01270166537925830792243875322356;
        %     0.88729833462074170213895740744192, 0.056350832689629148930521296279039;
        %     0.88729833462074170213895740744192, 0.1];
        % weights = [0.068464377671353521259689500766399, 0.10954300427416563401550320122624, 0.068464377671353521259689500766399, 0.061728395061728395061728395061728, 0.098765432098765432098765432098765, 0.061728395061728395061728395061728, 0.0086961161558069665494796751659123, 0.013913785849291147520001565851544, 0.0086961161558069665494796751659123];
    case "3d"
        points = [0.09273525035, 0.09273525035, 0.09273525035;
            0.72179424905, 0.09273525035, 0.09273525035;
            0.09273525035, 0.72179424905, 0.09273525035;
            0.09273525035, 0.09273525035, 0.72179424905;
            0.31088591930, 0.31088591930, 0.31088591930;
            0.06734224225, 0.31088591930, 0.31088591930;
            0.31088591930, 0.06734224225, 0.31088591930;
            0.31088591930, 0.31088591930, 0.06734224225;
            0.45449629590, 0.45449629590, 0.04550370415;
            0.45449629590, 0.04550370415, 0.45449629590;
            0.04550370415, 0.45449629590, 0.45449629590;
            0.45449629590, 0.04550370415, 0.04550370415;
            0.04550370415, 0.45449629590, 0.04550370415;
            0.04550370415, 0.04550370415, 0.45449629590];
        weights = [0.0122488405125, 0.0122488405125, 0.0122488405125, 0.0122488405125, 0.0187813209500, 0.0187813209500, 0.0187813209500, 0.0187813209500, 0.0070910034625, 0.0070910034625, 0.0070910034625, 0.0070910034625, 0.0070910034625, 0.0070910034625];

        % points = [0.25, 0.25, 0.25;
        %     0.319793627829630, 0.319793627829630, 0.319793627829630;
        %     0.040619116511110, 0.319793627829630, 0.319793627829630;
        %     0.319793627829630, 0.040619116511110, 0.319793627829630;
        %     0.319793627829630, 0.319793627829630, 0.040619116511110;
        %     0.091971078052723, 0.091971078052723, 0.091971078052723;
        %     0.724086765841831, 0.091971078052723, 0.091971078052723;
        %     0.091971078052723, 0.724086765841831, 0.091971078052723;
        %     0.091971078052723, 0.091971078052723, 0.724086765841831;
        %     0.443649167310371, 0.056350832689629, 0.056350832689629;
        %     0.056350832689629, 0.443649167310371, 0.056350832689629;
        %     0.056350832689629, 0.056350832689629, 0.443649167310371;
        %     0.056350832689629, 0.443649167310371, 0.443649167310371;
        %     0.443649167310371, 0.056350832689629, 0.443649167310371;
        %     0.443649167310371, 0.443649167310371, 0.056350832689629];
        % weights = [0.019753086419753, 0.011989513963170, 0.011989513963170, 0.011989513963170, 0.011989513963170, 0.011511367871045, 0.011511367871045, 0.011511367871045, 0.011511367871045, 0.008818342151675, 0.008818342151675, 0.008818342151675, 0.008818342151675, 0.008818342151675, 0.008818342151675];

        % points = [0.25, 0.25, 0.25;
        %     1/6, 1/6, 1/6;
        %     0.5, 1/6, 1/6;
        %     1/6, 0.5, 1/6;
        %     1/6, 1/6, 0.5];
        % weights = [-2/15, 3/40, 3/40, 3/40, 3/40];
    otherwise
        error("Invalid gauss-lengdre type.");
end
end